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Abstract. We study the mixing length concept in comparison to three-dimensional numerical calculations of 
convection with rotation. In a limited range, the velocity and temperature fluctuations are linearly proportional 
to the superadiabaticity, as predicted by the mixing length concept and in accordance with published results. 
The effects of rotation are investigated by varying the Coriolis number, Co = 2 fir, from zero to roughly ten, and 
by calculating models at different latitudes. We find that a decreases monotonically as a function of the Coriolis 
number. This can be explained by the decreased spatial scale of convection and the diminished efficiency of the 
convective energy transport, the latter of which leads to a large increase of the superadibaticity, 5 — V — V a d 
as function of Co. Applying a decreased mixing length parameter in a solar model yields very small differences 
in comparison to the standard model within the convection zone. The main difference is the reduction of the 
overshooting depth, and thus the depth of the convection zone, when a non-local version of the mixing length 
concept is used. Reduction of a by a factor of roughly 2.5 is sufficient to reconcile the difference between the 
model and helioseismic results. The numerical results indicate reduction of a by this order of magnitude. 

Key words, convection - hydrodynamics - stellar evolution 

1. Introduction to the fact that the temperature gradient there is nearly 

adiabatic, and the actual description of convection does 
The mixing length concept was originally developed to not make mudl difference . 

deal with terrestrial incompressible convection (e.g. Taylor 

HS3 PrandtllM. Later, this approach was adopted in 0ne wa ^ to test the validit y of the mixin S len S th con " 
the context of stellar convection (BiermannHHS Cowling 18 to P erf °™ numerical calculations of convection. 
mm . and the version introduced by Vitense <HM| see £2^P arlS0ns have J^ done e ' g ' by Chan & Sofia ( - 1987 ' 
also Bohm-VitenseESEl is still the most widely used de- ™> Kim et aL ™> and Porter & Woodward (2000). 
scription for convection in stellar structure models. The The mam result of these studies 18 that the relations de- 
main advantage of this model is its simplicity; only one rlved between varl0us mean thermodynamic and kme- 
fixed free parameter is needed which relates the correla- matic quantities under the basic mixing length assump- 
tion length to the local pressure scale height. In stellar tl0n are S ood approximations when the superadiabaticity, 
structure and evolution models the mixing length param- 5 = V ™ V °*> 18 sma11 enou S h ( Klm et aL conclude that 
eter, a, is needed in order to determine the temperature :t 18 S ^Sf nt to have 6 ~ °- 01) ' Intcrestln g!y, Chan & 
gradient in the convection zone. The value of a fixes the Sofia G2B3 And that the basic mixing length assumption, 
depth of the convection zone in the model, which, for the 1 = aH P> 1S realised m thcir model. 

solar case, is known from helioseismology (e.g. Monteiro The mixing length concept, and the aforementioned 

et al. 119941 Christensen-Dalsgaard et al. I1995f) . comparisons to numerical convection, all neglect the ef- 

Despite the crude nature of the mixing length concept, fects of rotation, which can be dynamically important in 

it can be useful in the bulk of the solar convection zone due many stars with convective envelopes. For example, the 

mixing length models of the solar convection zone (e.g. 
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approximate size of the convective elements to be of the 
order of the pressure scale height (H p 5 • 10 7 m), one 
can estimate the turnover time to be r w 10 6 s. The solar 
angular velocity is « 2.6- 10 _6 s _1 , with which one can 
estimate the Coriolis number, Co = 2 fir, to be of order 
ten near the bottom of the convection zone and of order 
unity in a large fraction of the convection zone (Fig.^). 

We present three-dimensional calculations of convec- 
tion in a rectangular domain that represents a small por- 
tion of a full star at some latitude. The latitude of the 
box can be chosen by imposing a suitably oriented ro- 
tation vector. In the present paper, our main aim is to 
parametrise the effects of rotation on the mixing length 
a. This is achieved by calculating several boxes with vary- 
ing rotational influence and latitude. One can think of the 
variation of the Coriolis number as the variation of the 
radial position in a convection zone. Thus, we obtain a 
radial profile for the mixing length a which we introduce 
into a solar model. 

In Sect.|3 of this paper the computational model is 
described, and in Sect.|3|a short summary of the mixing 
length concept is given. In Sects. 0] and 0we discuss the 
parameter ranges of the calculations and summarise the 
results. 

2. The model 

A description of the model setup can be found in Kapyla 
et al. I)2()()4al hereafter Paper I), and that of the numeri- 
cal method in Caunt & Korpi 1|2001[1 . The computational 
domain is a rectangular box, situated at a latitude O. The 
coordinates are chosen so that x, y, and z correspond to 
the south-north, west-east, and radially inward directions, 
respectively. The angular velocity as function of latitude 
is i~2 = f2(cos Qe x — sin 0e z ). 

The gas is assumed to obey the ideal gas equation 

p = pe(7-l), (1) 

where 7 = cp/cy = 5/3 is the ratio of the specific heats. 

To model the radiative losses near the surface we use 
a narrow cooling layer on the top of the convection zone, 
cooled with a term 
1 



^cool 



f(z)(e - e ) , 



(2) 



where i coo i is a cooling time, chosen to be short enough for 
the upper boundary to stay isothermal, f(z) a function 
which vanishes everywhere else but in the interval z < 
z < Z\, and eo = e(zo) the value of internal energy at the 
top of the box. 

We adopt periodic boundary conditions in the horizon- 
tal directions, and closed stress free boundaries at the top 
and bottom. The temperature is kept fixed at the top of 
the box and a constant heat flux is applied at the bottom 
du x du v „ . 

zo,z 3 ; (3) 



dz dz 







at z 



e(z ) 
de 
dz 



eo 



( 7 -l)(m 3 + l) 



(4) 
(5) 



where 7723 is the polytropic index in the lower overshoot 
layer. 

We obtain nondimcnsional quantities by setting 



Po 



cp 



1 



Thus, length is measured with respect to the depth of the 
unstable layer, d = z% — z\, and density in units of the ini- 
tial value at the bottom of the convectively unstable layer, 
po. Time is measured in units of the free fall time, yj d/g, 
velocity in units of \fdg~ , and entropy in terms of cp . The 
box has horizontal dimensions L x = L y = 4, and L z = 2 
in the vertical direction. The upper boundary is situated 
at zo = —0.15 and the upper boundary of the convectively 
unstable layer is at zi =0. The lower boundaries of the 
unstable layer and of the lower stably stratified region are 
at (z2,z s ) — (1,1.85), respectively. 

The dimensionless parameters controlling the calcu- 
lations are the Prandtl number Pr, and the Taylor and 
Rayleigh numbers, denoted by Ta and Ra, respectively. 
The relative importance of thermal diffusion against the 
kinetic one is measured by the Prandtl number 



Pr 



v 
Xo 



(6) 



where xo is the reference value of the thermal diffusivity, 
taken from the middle of the unstably stratified layer. In 
the present study, v is a constant. We set Pr = 0.4 which 
is a compromise between computational time (large diffu- 
sivity enforces a smaller time step) and realistic physics 
(in the Sun Pr < 1). 

Rotation is measured by the Taylor number 



Ta = 



/2fM 2 \ 2 
KIT) ' 



(7) 



A related dimensionless quantity is the Coriolis number, 
which is the inverse of the Rossby number, Co = 2 Or, 
where r = d/ut is the convective turnover time, and ut the 
rms-value of the fluctuating velocity determined through 
the calculation and averaged over the convectively unsta- 
ble layer and time. In our calculations Co varies by about 
two orders of magnitude from about 0.1 to roughly 14. 

Convection efficiency is measured by the Rayleigh 
number 



Ra 



d 4 gS 
Xo"H p 



(8) 



where S = V — V a d is the superadiabaticity, measured 
as the difference between the actual and adiabatic loga- 
rithmic temperature gradients, and H p the pressure scale 
height, both evaluated in the middle of the unstably strat- 
ified layer in the non-convecting hydrostatic reference so- 
lution. 

We define the Reynolds number as 



Re 



u t d 



(9) 



The initial stratification is polytropic, described by 
the indices mi, 7712, and 7713 for the three layers. We set 
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Table 1. Summary of the convection calculations. 
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Fig. 1. Logarithm of the Coriolis number, Co = 
2flQaHp/v, in the solar convection zone. 

mi = oo, to 2 = 1, and m 3 = 3, respectively, which means 
that the cooling layer is initially isothermal, and the strat- 
ification of the lower overshoot layer resembles that of the 
solar model of Stix l(2UU2jl . 

Initially the radiative flux, ,F ra d = «Ve, where k = 
JPX is the thermal conductivity, carries all of the energy 
through the domain. This constraint defines the thermal 
conductivities in each layer as 

* = ^*±± . (10) 

Kj rrij + 1 

In the calculations the vertical profile of k is kept con- 
stant, which with the boundary condition for the internal 
energy, Eq. J5J, assures that the heat flux into the domain 
is constant at all times. The initial state is perturbed with 
small scale velocity fluctuations of the order of 10~ 2 y/dg 
which are deposited in the convectively unstable layer. 

The code is parallelised using a message passing in- 
terface (MPI). The calculations were carried out on the 
IBM eServer Cluster 1600 supercomputer hosted by CSC 
Scientific Computing Ltd., in Espoo, Finland, and on the 
KABUL and BAGDAD Beowulf clusters with 16 and 34 
processors, respectively, at the Kiepenheuer-Institut fur 
Sonnenphysik, Freiburg, Germany. 

3. The mixing length concept 

The basic assumption in the mixing length concept is that 
the rising and descending convective elements preserve 
their identity for a length which is proportional to the 
local pressure scale height, i.e. 



/ = aH n 



(11) 



where H u 



l p = (pff ) " Using this assumption it is possible 
to derive equations for the mean convective velocity and 
temperature fluctuation (e.g. Stix[2002) 
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'T' 2 ) = ^(V-V ad )(T) , 



(13) 



where adiabatic variation of the convective eddies has been 
assumed, and the brackets denote horizontal and temporal 
averaging. 

The mixing length parameter enters the stellar evolu- 
tion models in the context of energy transport. Thus, we 
relate the fluxes of enthalpy and kinetic energy with the 
additional equations (e.g. Porter & Woodward 12000(1 



Fkin 



a c c p (p) 
"k(p)( 




(14) 
(15) 



(V- V 



ad l 



(12) 



where F c = (c p pT'u' z ) and Fki n = {\p u u'z) are the en- 
thalpy and kinetic energy fluxes, respectively. We use sub- 
scripts u, T, e, and k for the mixing length parameters 
obtained from the relations 1)12(1 to ((15(1 due to the fact 
that it is not a priori clear that all of them would yield 
the same value. 



4. Parameter ranges in the Sun and in the 
calculations 

4.1. Coriolis number in the Sun 

The Coriolis number varies strongly in the outer convec- 
tion zone of the Sun. Near the surface it is only of order 
10 -3 or smaller, but it exceeds 10 in the deeper region. 
Figure shows Co for a standard solar model, calculated 
according to Co = 2Q1/ v, where I = aH p is the mixing 
length, and v is the convection velocity that is obtained 
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Fig. 2. Correlations of u z with respect to the reference levels z rc { — (solid), 0.17 (dotted), 0.26 (dashed), 0.5 (dot- 
dashed), 0.77 (triple dot dashed), and 0.87 (long dashed) from the run s2CoO as function of the number of pressure 
scale heights, N p , (left) and the actual distance Az (right). 




-0.002 0.000 0.002 0.004 0.006 



Fig. 3. Relations between the vertical velocity and tem- 
perature fluctuations versus the superadiabatic gradient 
from the run s2Col-90. 



by means of the local mixing-length formalism. The co- 
efficient a = 1.66 is obtained through a calibration of a 
model of age 4.57 • 10 9 years to the present Sun. 

4.2. Summary of the calculations 

Table ^ gives a summary of the calculations. Bearing in 
mind the result presented in Fig. ^ we have made calcu- 
lations with four different Taylor numbers, which roughly 
correspond to Coriolis numbers 0, 0.1, 1, and 10. We 
refer to these runs with a prefix CoO, CoOl, Col, and 
ColO, respectively. If we interprete the Coriolis number 
as the depth dependence, the two first cases represent 
roughly the top, and the latter two the middle and the 
bottom of the solar convection zone. To investigate the 



latitudinal dependence we have made calculations at four 
latitudes, namely 9 = 0° (equator), —30°, —60°, and 
—90° (southern pole). In the standard setup the essentially 
step-function form of the thermal conductivity profile, see 
Eq. I|10|l , introduces a peak in the superadiabaticity at the 
bottom of the convection zone (not shown). We find this 
feature to persist in calculations where the resolution is 
doubled from the values listed in Table fusing runs from 
Paper I. 

Thus we introduce a setup in which the logarithmic 
temperature gradient in the lower overshoot layer in the 
initial state follows a profile 

V = V 3 + i{tanh[4(z m - z)] + 1}AV , (16) 

where V3 is the applied gradient at the bottom, and AV 
the difference between the gradient in the unstable layer 
and the applied gradient, and z m inflection point of the 
tanh-function, calculated so that V = V a( j at the interface 
between the stable and unstable layers. The thermal con- 
ductivity is determined so that initially radiative diffusion 
transports the total energy flux. We denote runs with this 
setup by a prefix s and a number which depends on the 
used resolution in Table ^ 

5. Results 

5.1. Basic mixing length assumption 

We test the basic mixing length assumption by calculating 
the correlation of the horizontally averaged vertical veloc- 
ity as function of depth as suggested by Chan & Sofia 

C\ 1 — ( u z rc f u z) (T7'\ 

^[U Zrc! ,U zi - (w 2^ )1/2(m 2 )1/2 , (U) 

where z re f is the reference level and z runs from zq to Z3. It 
is useful to plot this quantity as a function of the number 
of pressure scale heights 

N p = A\n(p), (18) 
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Fig. 4. The mixing length parameters calculated from Eqs. 

where the difference is calculated from the reference level. 
If the basic mixing length assumption, defined by Eq. ifTTl) , 
is correct the correlations plotted as functions of N p with 
respect to different reference levels should fall on top of 
each other. Figure |21 shows the correlations in the run 
s2CoO case as functions of N p and the distance measured 
in terms of d. The widths of the correlations measured as 
function of pressure scale heights decrease monotonically 
as a function of depth. We find that this width is more 
or less constant when measured in terms of d, the unit 
length. This result resembles the one by Robinson et al. 
(2003}, who present calculations of convection represent- 
ing the highly superadiabatic near-surface layers of the 
Sun. 

Although we do not find a scaling of the correlation 
width with the scale height reported by Chan & Sofia 
(1987, 1989), we do find that the velocity and tempera- 
ture fluctuations follow a more or less linear relation as 
function of the superadiabaticity, S, for a limited range of 
parameters (see Figs. 01 and 0}. This is realised in in the 
calculations where 8 decreases monotonically within the 
convectively unstable region due to the more smoothly 
varying thermal conductivity (Figs. and 0). 

The linear part in Figs. @ and (0J is clearly shorter 
than in the results of Chan & Sofia ifTWl [15531 HMfy 

which can be explained by the fact that whereas in our 
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calculations the convectively unstable region contains less 
than two pressure scale-heights, Chan & Sofia typically 
have five or more. Our failure to find any such linear rela- 
tion for the standard models is most probably due to an 
unrealistic step function like thermal conductivity profile 
which somewhat distorts the thermal structure. However, 
the qualitative trends seen as function of rotation do not 
depend on the choice of the thermal conductivity profile. 

5.2. Effect of rotation and dependence on latitude 

Figure 0] shows the mixing length parameters, as defined 
by Eqs. (|12(l to (|15|l . from the s2 runs with four Coriolis 
numbers at latitude —90°. These figures show that we do 
not find one universal parameter a that could describe all 
the quantities in Eqs. I|12|) to l|15[). Firstly, a e describing 
the convective energy transport, seems to have a rather 
well-defined value of about 1.5 for all cases bar the runs 
with largest rotation. 

Similarly, stays fairly constant for slow rotation, 
and shows a large reduction for the ColO case. On the 
other hand, cut and a u decrease rapidly as Co grows. The 
explanation to this behaviour is that the superadiabatic- 
ity, which appears in the denominator in the equations 
for «t and a u , increases rapidly as a function of rotation 
(see Fig. EJl . This growth of V — V ac i is linked to the de- 
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Fig. 5. Superadiabaticity V 
and rotation. 



V a d as function of depth 



crease of the convective energy transport as a function of 
rotation (see Paper I). The reduced convective flux forces 
the radiative diffusion to transport more energy, which 
can only be achieved by steepening the temperature gra- 
dient. However, due to the fact that the present numerical 
model deals with inefficient convection, where maximally 
about one third of the total energy flux is transported by 
convection, the rotational influence on the superadiabatic- 
ity is probably underestimated in comparison to the solar 
case where convection is expected to transport practically 
all energy. Preliminary numerical results with an efficient 
convection setup support this conjecture. 

Fig. shows the superadiabaticity in the overshoot 
layer from the s2 set of calculations. The increasing effect 
of rotation is seen in the larger magnitude of the superadi- 
abaticity and a significantly deeper superadiabatic layer. 
Furthermore, as the overshooting decreases as a function 
of rotation the transition at the unstable/stable interface 
becomes steeper due to the larger S in the convectively 
unstable layer. The quantitative change of 5 in the depth 
of the superadiabatic layer is probably exaggerated due 
to the much larger energy flux in the present calculations 
in comparison to the Sun. A similar issue concerns the 
transition from the overshoot to the radiative layer noted 
by Rempel (2004 ). According to that study, reducing the 
input flux sufficiently one should reach a regime where 
the overshoot resembles that obtained by non-local mixing 
length models, i.e. almost adiabatic overshoot and a sharp 
transition to the radiative gradient. However, this may re- 
quire the use of an anelastic code due to the timestep 
requirement. 

The main feature of the latitudinal dependence is that 
the superadiabaticity increases as function of latitude (see 
Fig. [3). Note that whereas the superadiabatic region ex- 
tends as the latitude increases, the overshooting depth de- 
creases in a similar fashion to Fig. |S| In an earlier study 
with denser latitude coverage (Paper I) some indications 
of a maximum of S at latitudes 9 « —60° ... — 75° was 



Fig. 6. Superadiabaticity V — V a d as a function of depth 
and rotation in the transition layer between the convec- 
tively unstable and stable regions. Note that the scale is 
different by roughly an order of magnitude in comparison 
to Fig. El 




Fig. 7. Superadiabaticity V — V a d as function of depth 
and latitude. 



observed. What this implies is that convection is more 
efficient at low latitudes and that due to this one would 
expect overshooting to peak there. This is indeed observed 
for slow and moderate rotation, i.e. Co < 4, see Fig. 9 of 
Paper I. This could mean that the overshooting extends 
deeper at equatorial regions, which can lead to a prolate 
shape of the tachocline as indicated by helioseismology 
(Basu & Antia |2"001II . However, one must bear in mind 
that the helioseismology results refer to the shear layer 
whose connection with overshooting is not yet known. The 
thermal stratification of the Sun does not vary observably 
as a function of latitude as in the numerical models which 
is most likely due to the much too vigorous convection in 
the calculation. Furthermore, we note that the trend in 
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Fig. 8. Mixing length parameters a;, defined via Eqs. Ijl2(l to l|15(l . as function of Ra and resolution for the runs siCoO 
from Tabled 



the overshooting depth as function of latitude seems to 
reverse for more rapid rotation (Paper I). 

We feel that further study, beyond of the scope the 
present paper, of the details of the overshooting as func- 
tions of rotation and input energy flux is needed in order 
to determine whether it is possible to approach the mix- 
ing length regime with 3D calculations and to fully sub- 
stantiate the effects of rotation which are virtually always 
neglected in overshooting models. 

5.3. Effects of resolution 

In order to study the dependence of the mixing length re- 
lations on resolution, we have made a set of runs where 
we vary the Rayleigh number from 3.1 • 10 5 to 4.1 • 10 6 and 
the number of gridpoints from 64 3 to 256 3 , named slCoO 
to s4CoO in Table d Fig. |H1 shows the mixing length pa- 
rameters from the runs slCoO to s4CoO. The main feature 
of the results is that the qualitative character of the a- 
parameters remains unaffected when resolution changes. 
The relatively large changes in the parameters and a u 
are mostly explained by the smaller superadiabaticity as 
resolution increases. 

The largest resolution calculation, however, is not fully 
comparable to the lower resolution ones due to the fact 



that the calculation did not reach a thermally saturated 
state. This has to do with the fact that in order to avoid 
the long thermal relaxation phase we start the calculation 
with a stratification that is expected to be close to the final 
thermally relaxed state. This procedure has to be carried 
out by trial and error for each resolution and for the 256 3 
run the initial guess of the stratification was too superadi- 
abatic, resulting in too efficient convection. However, we 
think that even though this run is not quantitatively fully 
comparable to the others, it still illustrates the fact that 
the qualitative results remain unchanged as a function of 
resolution. 



5.4. Application to a solar model 

As the Coriolis number in the lower part of the convection 
zone is large, the convective efficiency should be substan- 
tially reduced according to the 3D numerical calculations. 
In the one-dimensional standard solar model we simulate 
this by a reduced mixing length. However, there are vir- 
tually no consequences to the solar structure within the 
convection zone. This is because the temperature gradi- 
ent, as determined by mixing- length theory, closely fol- 
lows the adiabatic gradient, irrespective of the details of 
the convection formalism. 
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Fig. 9. Temperature gradients V = d\nT/d\nP near the 
base of the convection zone. The actual gradient (heavy 
line) is identical to the radiative gradient in the core, and 
closely follows the adiabatic in the convection zone and in 
the region of overshooting. A local model with no over- 
shooting (lower line); non-local models with mean path 
1/2 (above), as assumed in local theory, and with mean 
path 1/5 (middle). 

On the other hand, the overshooting below the base of 
the convection zone may be severely reduced by a reduced 
convective efficiency This can be illustrated by solar mod- 
els that employ a non-local version of the mixing-length 
formalism (Skaley & Stix lTMTl Stix 2002 ). In those models 
the temperature excess AT and the convection velocity v 
are determined as integrals over the distance a convection 
parcel travels in the vertical direction. If this distance is re- 
duced, the parcel acquires less kinetic energy, so that there 
is less overshooting. In Fig. [5] this is shown for a standard 
model calculated with local mixing-length theory and two 
non-local models. We find that the depth of overshooting 
decreases approximately in proportion to the decrease of 
the integration path. If the average parcel path is reduced 
by a factor 2.5, which is consistent with the numerical 
results, overshooting extends only over w 2900 km. 

Christensen-Dalsgaard et al. Ill995fl determined an up- 
per limit of rj O.li/p for overshooting with a sharp transi- 
tion from the adiabatic to the radiative regime, such as we 
obtain with the non-local mixing length models. Beyond 
this limit the discontinuity of the temperature gradient 
would leave a measurable signature in the spacing of p- 
modc frequencies. Thus, the rotational quenching of the 
overshooting would appear to be welcome. On the other 
hand, the total depth of the convection zone may become 
too small in this way. The problem is aggravated as a 
recent redetermination of the abundance of oxygen and 
other heavy elements (Asplund et al. 2004| lowers the 



opacity and, therefore, also leads to a shallower convection 
zone. Bahcall et al. (PUS) find r/r Q = 0.726 for the base 
of the convection zone of their model BP04+ which in- 
corporates the new lower opacity as well as best available 
updates of all other input parameters. This means that 
the convection zone would be too shallow by w 9000 km 
or sa 0.2H p , as compared to the helioseismological result 
(r/r Q = 0.713±0.001; Christensen-Dalsgaard et al. HMTl 
Basu & Antia [Tnn3)- 

Presently is is not clear how this conflict will be re- 
solved. New opacity calculations (Seaton & Badncll 2004) 
yield an increase of 5 % compared to the latest OPAL 
values; but this is only about half the increase needed: 
Bahcall ct al. ( 2005) obtain a sufficiently deep convection 
zone if they increase the opacity by s» 10 % in the depth 
range 0.4r Q to 0.7r Q . We suggest that the solution is a 
combination of modest overshooting and a modest opac- 
ity increase. Overshooting may also account for the whole 
difference if it avoids the sharp transition that occurs in 
our calculation. Rempel (2004} has proposed a model in 
which the mixing between downflows and upflows is cru- 
cial for the total depth of the overshooting, and in which 
an ensemble of downflows with different strength produces 
a smooth transition. 
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